Fast Fourier Transform Notes 
18.310, Fall 2005, Prof. Peter Shor 



1 Introduction: Fourier Series 

Early in the Nineteenth century, Fourier studied sound and oscillatory motion and conceived 
of the idea of representing periodic functions by their coefficients in an expansion as a sum of 
sines and cosines rather than their values. He noticed, for example, that you can represent 
the shape of a vibrating string of length L, fixed at its ends, as 

oo 

y(^) = '^(^k sm{Trkx/L) 

x=l 

The coefficients, a^, contain important and useful information about the quality of the 
sound that the string produces, that is not easily accessible from the ordinary y = f{x) 
description of the shape of the string. 

This kind of representation is called a Fourier Series, and there is a tremendous amount 
of mathematical lore about properties of such series and for what classes of functions they 
can be shown to exist. One particularly useful fact about them is the orthogonality property 
of sines: 

/ sm{'!rkx / L) sm{Tr jx/L)dx = Sj}- — . 
Jx=o ' 2 

Here Sj^k is the Kronecker delta function, which is if j 7^ A; and 1 if j = k. The integral 

above, then, is unless j = k. in which case it is L/2. To see this, you can write the 

product of these sines as a constant multiple of the difference between cos(7r(fc + j)x/L) 

and cos(7r(/c — j)x/L), and realize that unless j = ±k, each of these cosines integrates to 

over this range. 

By multiplying the expression for y{x) above by sin(7rja;/L), and integrating the result 
from to L, by the orthogonality property everything cancels except the sm{Trjx/L) term, 
and we get the expression 

2 

dj = — / f{x)sm[Trjx/L)dx 
L Jx=o 

Now, the above sum of sines is a very useful way to represent a function which is at 
both endpoints. If we are trying to represent a function on the real line which is periodic 
with period L, it is not quite as useful. This is because for a periodic function, we need 
/(O) = f(L) and /'(O) = f'{L). For the sum of sines above, the terms with odd k such as 
siuTTX are not themselves periodic with period L. For periodic functions, a better Fourier 
expansion is 

00 00 
y{x) = ao + ^ aj cos{2'it jx/L) + ^ 6^ sin(27rA;x/L). 

j=i k=i 
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It is fairly easy to rewrite this as a sum of exponentials, using the identities 
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This results in the expression (with a different set of coefficients aj) 



oo 



(1) 



j=-oo 



The orthogonality relations are now 



Jx=0 



.2-n-ij/L^2Trik/L _ ^ 



This means that we now can recover the aj coefficient from y by performing the integral 



2 The Fourier transform 

Given a function f(x) defined for all real x, we can give an alternative representation to it 
as an integral rather than as an infinite series, as follows 



Here g{x) is called the Fourier transform of f{x), and f{x) is the inverse Fourier transform 
of g{x). This is a very important tool used in physics. 

One reason for this is that exponential functions e*'^^, which / is written as an integral 
sum of, are eigenfunctions of the derivatives. That is, the derivative, acting on an expo- 
nential, merely multiplies the exponential by ik. This makes the Fourier transform a useful 
tool in dealing with differential equations. 

Another example of the Fourier transform's applications can be found in quantum me- 
chanics. We can represent the state of a particle in a physical system as a wave function 
4>{x), and the probability that the particle in this state has a position lying between x and 
X + dx is \(f){x)\^dx. 

The same state can also be represented by its wave function in momentum space, and 
that wave function of the variable p is a constant multiple of the Fourier transform of (f){x): 




(2) 
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We can derive a formula for computing the Fourier transform in much of the same way 
we can compute Fourier series. The resulting formula is 

where the integration is over all real values of x. 

The Fourier transform can be obtained by taking the Fourier series and letting L go to 
DO. This turns both the function and its Fourier series into functions defined over the real 
line. The finite Fourier transform arises by turning these both into a finite sequence, as 
shown in the next section. 

3 The Finite Fourier Transform 

Suppose that we have a function from some real-life application which we want to find the 
Fourier series of. In practice, we're not going to know the value of the function on every 
point between and L, but just on some finite number of points. Let's assume that we 
have the function at n equally spaced points, and do the best that we can. This gives us 
the finite Fourier transform. 

We have the function y{x) on points jL/n, for j = 0, 1, . . . , n — 1. We would like to 
represent the function 

k 

but, of course, we only have knowledge of y at values of x = jL/n. If we plug in x = jL/n 
for the value of x, we get 

k 

There is one more thing to notice here. Namely, that adding n to k doesn't change any 
of the values of e^'^*-'*^/" because e^'^* = 1. Thus, if we just have n equally spaced points, 
we only need the first n terms from the Fourier series to perfectly match the values of y at 
our points. This makes sense — if we start with n complex numbers yj, we end up with n 
complex numbers a^, so we keep the same number of degrees of freedom. The here are 
the finite Fourier transform of the yj. 

How do we compute the a^, given the yj. It's not hard to see that it works in essentially 
the same way that it did for the complex Fourier series we talked about earlier, only we 
have to replace an integral with a sum. Thus, we get 

^k = lT.yje-'^''''/\ (4) 
j 

Here, again the are the finite Fourier transform of the yj. The inverse Fourier 
transform given above (Eq. 3) is the same as the finite Fourier transform, except we have 
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put a — sign on the exponent, and left out the factor of 1/n. In fact, you will sometimes 
see the factor of 1/n distributed equally, with a 1/ -y/n on both the forwards and the inverse 
Fourier transforms. 

How do we prove that the formulas (4) and (3) are inverse transforms of each other? 
The proof works the same way as it does for the Fourier series, and in fact this formula can 
be derived from (2) and (1). The orthogonality relations turn into the sum 

n— 1 

i=o 

I'll let you figure out the rest of it. 

4 Computing the finite Fourier transform 

It's easy to compute the finite Fourier transform or its inverse if you don't mind using 
O(n^) computational steps. The formulas (3) and (4) above both involve a sum of n terms 
for each of n coefficients. However, there is a beautiful way of computing the finite Fourier 
transform in only 0{nlogn) steps. 

One way to understand this algorithm is to realize that computing a finite Fourier 
transform is equivalent to plugging into a degree n — 1 polynomial at all the n roots of 
unity, e^TTjfc/n^ ^^j. g < A: < n — 1. The Fourier transform and its inverse are essentially the 
same for this part, the only difference being which n-th root of unity you use. So, let's be 
consistent with Prof. KIcitman's notes and do the inverse Fourier transform. 

Suppose we know the values of and we want to compute the yj using the inverse 
Fourier transform, Eq. (3). Let the polynomial p{x) be 

n-l 
fe=0 

Now, let z = e^'^*/". Then, it is easy to check that we have 

This shows we can express the problem of the inverse Fourier transform as evaluating the 
polynomial p at the n-th roots of unity. 

What we will show is that if n is even, say n = 2s, it will be possible to find two degree 
s — 1 polynomials, Peven and Podd; such that we get all n of the values yj for < j < n — 1 
by plugging in the s-th roots of unity into Peven and Podd- The even powers of z will appear 
in Peven ) Sbud the odd powers of z will appear in Podd- If n is a multiple of 4, we can then 
repeat this step for each of Peven and Podd, so we now have our n values of yj appearing as 



4 



the values of four polynomials of degree n/4 — 1, when we plug the ^th units of unity, i.e., 
the powers of z^, into all of them. If n is a power of 2, we can continue in the same way, 
and eventually reduce the problem to evaluating n polynomials of degree 0. But it's really 
easy to evaluate a polynomial of degree 0: the evaluation is the polynomial itself, which 
only has a constant term. So at this point we will be done. 

The next question we address is: how do we find these two polynomials Pevcn and Podd? 
We will do the case of pcvcn first. Let us consider an even power of say z"^^ . We will look 
at the j-th term and the (j + s)-th term. These are 

ajz'^^^ and aj+.z^^^+^fe^ 



But since = = 1, we have 



^2kj+2ks _ ^2kj 



. Thus, we can combine these terms into a new term in the polynomial Pevenj 

with coeffi- 
cients 

= + %-l-s 

If we let 

s-l 

Peven — ^ ^ ^i'^"' 
j=0 

we find that 

=Peven(^^'')• 
Now, let US do the case of the odd powers. Suppose we are evaluating p at an odd power 
of z, say z'^^^^. Again, let's consider the contribution from the j'-th and the (j + s)-th terms 
together. This contribution is 

Here we find that 2^^'^+^)* = —1. Why? Again, z^*^* = 1, and z^ is a square root of 1, 
because when we square it we get z'^s = z^ = 1. Since ^; is a primitive n-th root of 1, z^ is 
not 1, so it must be —1. We now have 



Setting the j-th coefficient of Podd to 

^3 — ~ ^j+s 



and letting 

s-l 

Podd = X] 

j=0 
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we see that 

So now, we can show how the Fast Fourier transform is done. Let's take n = 2*. Now, 
consider an n x t table, as we might make in a spreadsheet. Let's put in our top row 
the numbers oq through a^-i- In the next row, we can, in the first ^ places, put in the 
coefficients of Peven, and then in the next ^ places, put in the coefficients of Podd- In the next 
row, we repeat the process, to get four polynomials, each of degree f — 1- After we have 
evaluated the second row, we treat each of Pcvcn and Podd separately, so that nothing in the 
first ^ columns subsequently affects anything in the last § columns. In the third row, we 
will have in the first | places the coefficients of Peven,even! which give us the value of p{z^^) 
when we evaluate Peven even{z^^)- Then in the next f places, we put in the coefficients of 



4k-] 



Peven.odd- This polynomial will give the value of p{z^^^'^) when we evaluate Pcvcn,oddl-2' 
The third ^ places will contain the coefficients of Podd,even, which gives us the values of 
p{z^^^^). The last | places will be occupied by the coefficients of Podd,odd) which gives the 
values of p{z'^^^^). Prom now on, we treat each of these four blocks of ^ columns separately. 
And so on. 

There are two remaining steps we must remember to carry out. The first is that the 
values of p{z^) come out in the last row in a funny order. We have to reshuffle them so 
that they are in the right order. I will do the example of n = 8. Recall that in the second 
row, the polynomial po, giving odd powers of followed pe-, giving even powers of z. In the 
third row, first we get the polynomial giving z^^ , then z^*^"*"^, then z^^~^^, then z'^^'^^. So 
in the fourth row (which is the last row for n = 8), we get the values of p{z'^) in the order 
indicated below. 



1 


2 3 


4 5 


6 7 


coefficients of p 


Pe{z'^'^) = p{z'^'') 




p,,,{z''^)=p{z^>') 


PeAz^'^)=piz^'^'') 








P{z') p{z^) 


P{z') p{z^) 


p{z-^) p{z^) 



You can figure out where each entry is supposed to go is by looking at the numbers in 
binary, and turning the bits around. For example, the entry in the 6 column is p{z^). You 
can figure this out by expressing 6 in binary: 110. You then read this binary number from 
right to left, to get Oil, which is 3. Thus, the entry in the 6 column is p{z^). The reason 
this works is that in the procedure we used, putting in the even powers of z first, and then 
the odd powers of z, we were essentially sorting the powers of z by the I's bit. The next 
row ends up sorting them by the 2's bit, and the next row the 4's bit, and so forth. If we 
had sorted starting with the leftmost bit rather than the rightmost, this would have put 
the powers in numerical order. So, by reversing the bits and sorting, we get the right order. 

The other thing we have to do is to remember to divide by n if it is necessary. We only 
need do this for the Fourier transform, and not the inverse Fourier transform. 
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5 Multiplication and Convolution 



Let's go back to the complex Fourier series. Suppose we have two functions / and g. 
Suppose we have the Fourier series for these two functions, that is, 

oo 

f{x)= Yl a,e2-i-/^ 

j=-oo 

and 

oo 

g{x)= Y ^fce^-'"/'' 

fc=— oo 

How do we find the Fourier series of the sum of these two functions? It's easy. We take the 
sum of the Fourier coefficients. 



fix)+gix)= Y iak + bk)e^^'''^/'' 

fc=— oo 

How about the product? This requires some calculation. 



oo 



f{x)9{x)= Y a^he^-'^^+^>/'^ 
j,k=-oo 

Now, what is the e^'^*'^/^ term of this Fourier series? We get a contribution to I exactly 
when j + k = I. So, 

oo / oo \ 

f{x)9{x) = E E ^A-j e^"''-/^ 

l=-oo \j=-oo / 

This sequence, q = Y^j ajbi^j, is known as the convolution of sequences aj and b^. We 
recently saw how convolutions come up in generating functions. If you multiply two gener- 
ating functions together, you take the convolution of their respective sequences. The same 
thing is true in the Fourier transform, in Fourier series and the finite Fourier transform: 
taking the Fourier transform turns pointwise multiplication turns into convolution, and vice 
versa. 

I didn't do this in class, but I think it might be worth going through in these notes. 
Let's check that the convolution of functions turns into multiplication of the Fourier series. 
Remember that f{x) and g{x) are both periodic with period L. Their convolution is 

h{y) = [ f{x)g{y - x)dx 

Jo 
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k 

where we only get the terms with j = k because Jq e^'^^^^~^^^/^ is if j 7^ k. 

Wait a minute! Where did that L come from? We didn't get it when we multiphed / 
and g. I don't have any good intuitive explanation for it right now, but if you look at the 
equations, you can see it's really there. Remember this extra L. We'll see it again in a 
couple of paragraphs. 

Now, let's work out the details of the fact that pointwise multiplication turns into 
convolution for finite Fourier series. Suppose we have two sequences and their finite Fourier 
transforms, i.e.. 



^2mjk/n 
3 

9k = ^6,e2-^'=/- 



What do we get when we take the inverse Fourier transform of the pointwise product /jt^fe? 
Let 



fk9k = X! 



3 

Then taking the Fourier transform, we get 



Cl = -Y.^-'''''"''fk9k 



" k 



= 1 ^ Q-2'Kilk/n ^ ^,^2mjk/n ^ ^^^^27Tij'k/n 



k J 

1 



2mk{j+j'-l)/n 

n 



3 3' k 

= ^(^3^1-3 

3 

where the last equality holds because the sum over fc is unless I = j + j'. Thus, pointwise 
multiplication turns into convolutions for the finite Fourier transform as well. 

We can use this fact that pointwise multiplication turns into convolution to multiply 
polynomials efficiently. Suppose we have two degree d polynomials, and we want to multiply 
them. This corresponds to convolution of the two series that make up the coefficients of 
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the polynomials. If we do this the obvious way, it takes 0{d^) steps. However, if we use the 
Fourier transform, multiply them pointwise, and transform back, we use 0{dlogd) steps for 
the Fourier transforms and 0(d) steps for the multiplication. This gives 0{dlogd) total, 
a great savings. We must choose the n for the Fourier series carefully. If we multiply two 
degree d polynomials, the resulting polynomial has degree 2d, or 2d + 1 terms. We must 
choose n > 2(i+ 1, because we need to have room in our sequence ag, ai, ... dn-i for all the 
coefficients of the polynomial; if we choose n too small, the convolution will "wrap around" 
and we'll end up adding the last terms of our polynomial to earlier terms. 

One can ask the question: when you multiply polynomials by taking the Fourier trans- 
forms and multiplying pointwise, how many times do you have to divide by n? At first 
sight, it seems like the answer should be inconsistent. If wc first take the inverse Fourier 
transform, Eq. (3), of each of the polynomials, multiply then, and then take the Fourier 
transform, we end up dividing by n when we take the forwards Fourier transform, so we 
divide by n once. On the other hand, if we take the forwards Fourier transform of the poly- 
nomials first, we have to divide each of these polynomials by n. Then, when we multiply 
them together, we have divided each of these by n, so we've ended up dividing by total. 
Sot it seems like, when we take the inverse Fourier transform, we should end up with an 
extra division by n than we did if we took the inverse Fourier transform first. Both of these 
methods can't be right, since the two answers differ by a factor of n. What's going on? 

The correct answer, as you should be able to figure out if you sit down and write 
everything out, is that you only divide by n once. This is related to the extra L term we 
mentioned above, which turns into an n when you turn Fourier series into finite Fourier 
transforms, and cancels out one of the extra factors of n. 

6 Fourier transforms modulo p and fast integer multiplica- 
tion 

So far, we've been doing finite Fourier transforms over the complex numbers. We can 
actually generalize them to work over any field with a primitive n-th root of unity. Suppose 
z is our n-th root of unity. The main fact we need to have the finite Fourier transforms 
work is that if z ^ and = 1, then 

n-1 

A;=0 

This holds for any field with a primitive n-th root of unity. The transforms work the same 
way as in Eqs. (3) and (4). The inverse Fourier transform is 

n-1 
k=0 
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and the forward Fourier transform is 

2 n—l 

If we take a prime p, then the field of integers mod p has an n-th root of unity if 
p = mn + 1 for some integer m. In this case, we can take the Fourier transform over 
the integers mod p. Thus, 17 has 16th roots of unity, one of which can be seen to be 
3. So if we use 2; = 3 in our fast Fourier transform algorithm, and take all arithmetic 
modulo 17, we get a finite Fourier transform. We can use this for multiplying polynomials. 
Suppose we have two degree d polynomials, each of which has integer coefficients of size at 
most B. The largest possible coefficient in the product is (5 — l)^((i + 1). If we want to 
distinguish between positive and negative coefficients of this size, we need to make sure that 
p > 4(i? — l)^((i+l). We also need to choose n > 2d+l, so as to have at least as many terms 
as there arc coefficients in the product. We can then use the Fast Fourier transform (mod p) 
to multiply these polynomials, with only 0{dlogd) operations (additions, multiplications, 
taking remainders modulus p), where we would have needed d"^ originally. 

Now, suppose you want to multiply two very large integers. Our regular representation 
of these integers is as Y.^'^k^^^^ where d^ are the digits. We can replace this by Ylik'^k^^ 
to turn it into a polynomial, then multiply the two polynomials using the fast Fourier 
transform. 

How many steps does this take? To make things easier, let's assume that our large 
integers are given in binary, and that we use a base B which is a power of 2. Let's assume 
the large integers have m bits each and that we use a base B (e.g., 10 in the decimal system, 
2 in binary) that has b bits. We then have our number broken up into m/b "digits" of b bits 
each. How large docs our prime have to be? It has to be larger than the largest possible 
coefficient in the product of our two polynomials. This coefficient comes from the sum of 
at most m/b + 1 terms, each of which has size at most (2^ — 1)^ < 2^^ This means that 
we're safe if we take p at least 

(y + 1)2^^ 

or taking logs, p must have around 26 + log2 y bits. 

Rather than optimizing this perfectly, let's set b = log2 m] this is simpler and will give 
us the right asymptotic growth rate. We thus get that p has around 31og2m bits. We 
then set n = 2^ + 1, so that our finite Fourier transform involves 0(n log n) = 0{m) 
operations, each of which may be an operation on a (31og2m)-bit number. If we use 
longhand multiplication and division (taking 0{b^) time) to do these operations, we get an 
0(m log^ m)-time algorithm. 

There's no reason that we need to stop there. We could always use recursion and 
perform these operations on the 36-bit numbers using fast integer multiplication as well. 
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If we use two levels of recursion, we get an (^(log n(log log n)^) time algorithm. If we use 
three levels of recursion, we get an O (log n (log log n) (log log log n)^ time algorithm, and so 
forth. 

It turns out, although I won't go into the details, that if you work hard enough you can 
get a 0(n log n log log n) time algorithm. The details can be found in Aho, Hopcroft and 
Ullman's book "Design and Analysis of Computer Algorithms." This algorithm uses the 
Chinese remainder theorem to get a recursion that produces a running time just a little bit 
more efficient than what you get using recursion with the method above. 

7 Details of making a spreadsheet 

It seems at first like the FFT is perfectly suited for a spreadsheet. It can be illustrated 
beautifully by a two-dimensional diagram. However, once you start looking at the details, 
it gets fairly tricky. The diagram illustrating the FFT is: 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 




FFT Figure 



Here, the outlined boxes in the top row contain our input. In every black box we put 
On + o,i^s where 04 and a^-i-^ are the entries in the two boxes in the row connected to this 
black box by black lines. In each yellow box, we put {a^ — a^^s)z^, where now a, and a,^^ 
are the boxes in the row above connected to this yellow box by brown lines. One thing to 
remember is that z (as well as s) changes from row to row. In calculating the boxes in the 
second row, we use our original z. In the third row, we use z"^. (So in terms of our original 
z, we are really multiplying by 2;^* here and not z^.) In the fourth row, we use z'^. In the 
fifth row, it turns out we don't need to use z at all, since in all these yellow boxes we have 
i = 0. 

One tricky thing about this is that we have to start counting i from the first box of the 
contiguous row of yellow boxes it's in, and not from the column. So, for example, in row 3 
you need to get i by taking the column number mod 4 and not just the column number. 
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After row 5, we're done with all the addition steps, and have only one or two steps left. 
In row 6, we shuffle all the numbers into the red boxes according to the reversal-of-binary- 
digits permutation. If you're constructing the spreadsheet according to the instructions in 
Prof. Kleitman's lecture notes, you don't actually need row 6, because he folds this shuffle 
into rows 2 through 5. However, you haven't really gotten rid of any complexity, because 
he has made constructing rows 2 through 5 arc a little more complicated. 

Finally, depending on whether we're doing the FFT or the inverse FFT, we may need 
to make a row 7 where we divide row 6 by n, which in this case is 16. 

If we're doing an FFT mod p for some prime all of these arithmetic operations are 
going to have to be done mod p, and in the last step, we will have to multiply by n~^, 
which is the integer satisfying 

n ■ = 1 mod p 

This inverse can be found by Euclid's gcd algorithm. 

There are a number of issues that can arise when you construct the spreadsheet for the 
homework. The first is the issue of copying formulas in spreadsheets. Suppose you want 
to give a spreadsheet for the FFT in which you can easily change the prime p and/or the 
root of unity z. If you put z oi p into the formulas, then when you copy the spreadsheet 
you'll have to change them all by hand. Furthermore, if you get them by using an absolute 
reference like $A$3 (say to a cell containing p or z) then when you copy this cell, the 
reference in the copy will still points to A3. If you want to change the values of p and z, 
you will need to replace all these references. 

Of course, you can argue that this is what the find-and-replace tool in the spreadsheet 
editor is intended for. But Prof. Kleitman has also figured out how to construct a spread- 
sheet cleverly so that we don't have to use this tool. Suppose we start our spreadsheet with 
the three rows 





A 


B 


C 


D 


E 


F 


G 


H 


I 


J 


K 


L 


1 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 ... 


2 


p 


p 


p 


p 


p 


p 


p 


p 


p 


p 


p 


p ... 


3 


z 


z 


z 


z 


z 


z 


z 


z 


z 


z 


z 


z ... 



Then if we want to make a reference to p, we can use a reference like C$2. This uses 
the value in the same column, but the second row. Now, we can copy our spreadsheet 
horizontally. This changes C to column farther to the right, and we can put a new prime 
p in the second row of our copy. As long as we don't change the rows our spreadsheet 
occupies, but just move it horizontally, we're fine. 

We now need to calculate the values of hi using the equations 

bi = CLi + di+s 



12 



The first case (bi) is easy. The second case (6,) is trickier, because of the term. The 
obvious way is to calculate it is to use a statement like 

= MOD( (J5 - N5) * z^^*'\ p). 

to calculate N6 in row 6, where s = 4 and we multiply by {z'^y. There are two problems 
with something like this. The first is that, if row 2 is all p's and row 3 is all z's, we can get 
p and z by N$2 and N$3, but where do we get the i from? In row 3, we have two stretches 
of four yellow cells, and we want i to start counting from each time. If we just reference i 
up in the first stretch, we'll have to change the formula if we copy it to the second stretch. 
This can be solved by using MOD(z, 4) in the formula. 

The second problem is somewhat more serious. This is that using the spreadsheet's 
arithmetic to compute z^* is very likely to give you an integer overflow, if z or i is large. 
What we could do is to create an extra row that contains powers of z mod p, and use the 
OFFSET function to find the right cell in it. This works fine, but it's kind of complicated 
and messy, and it's easy for mistakes to sneak in. 

Here's what I think might be the easiest solution: create another (n x log2 n) array in 
your spreadsheet that contains the numbers that you need to multiply by — a^+s in 
the h cells. For n = 16, the numbers you need to multiply them by look like 









- 1 


z 


z^ 


z^ 


z^ 


z^ 


z6 


z^ 




1 


z^ 


^6 - 








1 




z^ 


z^ 


- - 1 z^ 




1 


z^ - 




1 


z^ 






1 


z^ 


-1-1 


- 1 




1 - 


1 




1 




1 




1 



where all of the powers are performed mod p. Here — means that we have b instead of b in 
the corresponding cell, so there's no term needed for that cell. 

How can we constrct this second array? Since we'll never actually be looking at the — 
cells, we can fill them in anyway we want, and it's easiest to fill them in with the same 
sequence we find in b cells. 

For example, for n = 16, if p = 193, and z = 3 is our 16th root of 1, we get 



1 


3 


9 


27 


81 


50 


150 


64 


1 


3 


9 


27 


81 


50 


150 


64 


1 


9 


81 


50 


1 


9 


81 


50 


1 


9 


81 


50 


1 


9 


81 


50 


1 


81 


1 


81 


1 


81 


1 


81 


1 


81 


1 


81 


1 


81 


1 


81 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 
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How do we create this table? It's easy. We get the first eight cells in the top row by sticking 
in a 1 on the left, and then repeatedly multiplying by 3 mod 193. We can get the rest of 
the cells by copying the cell 8 positions to the left. In the second row, we can get the first 
four cells by squaring the cells above them mod 193, and the rest by looking 4 positions to 
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the left. In the third row, we can get the first two cells by squaring the celle above them, 
and the rest by looking 2 positions to the left. Larger FFT's will work similarly. 

Unless you want to use a search-and-replace function, you should remember to get 3 
and 193 by referencing cells that contain them, rather than coding them into the formulas. 

Finally, we need to shuffle all the entries around, by copying the entry in each cell to the 
column obtained by reversing the digits of the column's number in binary. For example, 
if n = 16, you would move the entry in the 5th (= 0101) position to the 10th (= 1010) 
position. I didn't find it too painful to construct this shuffle by hand. I also don't see 
any easy way of making the spreadsheet do it for you, but I wouldn't be too surprised if 
someone thinks of a clever way. 

You also need to remember to divide by n. You do this by multiplying by n^^ (mod p). 
Calculating n"^ can be done using the Euclidean algorithm, which you should have gone 
over earlier in 18.310. 

There are a lot of ways to make errors in this process. One thing you can do to make 
it easier is to start out by making a spreadsheet which uses z = 3 and p = 17 or 193, and 
change the numbers to be bigger later. You can also put in test sequences where you can 
easily figure out what the FFT should be, such as putting in one 1 and the rest O's. I find 
it useful to color the cells using the b formula one color and those using the h formula a 
different color, so that I could visually see if I stuck these formulas in the wrong columns 
by mistake. 

Fve explained how to do the FFT on a spreadsheet using two different types of cells 
in each row: one type that computes h and another that computes h. If you use slightly 
more complicated formulas, you might be able to tisc jiist two formulas total for all of these 
rows, again, one for h and one for 6, but where these two types don't need to be customized 
for each rows. Here, I don't see any way to do it without using use the OFFSET function 
(which generally makes things more complicated) , and you also need some way to figure out 
which row you're in (which also makes things more complicated, although not by much), 
so unless you can figure out some clever way of doing it without the OFFSET function, it's 
probably not worth the extra trouble. 
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